#include <iostream>
#include <fstream>
#include <vector>
#include <math.h>

#include "PGPLOT_Renderer.h"
#include "DataContainer.h"

#include "cpgplot.h"
#include "vis_utils.h"

void zphi2vec( double z, double phi, vector<double>& vec)

{ 
	double x,y,sz;

	x = cos(phi);
	y = sin(phi);
	sz = sqrt((1. - z)*(1. + z));

	vec[0] = sz*x;
	vec[1] = sz*y;
	vec[2] = z ; 

}


PGPLOT_Renderer::PGPLOT_Renderer() {

	cout << "PGPLOT_Renderer::PGPLOT_Renderer Default Constructor" << endl;
	
	this->device = NULL;
	this->verbose = false;
	xmin=-43200;
	xmax=43200;
	ymin=-324000.;
	ymax=324000.;

	hi_limit = -1;
	low_limit = -1;


	title = "";
	xlabel = "Right ascension";
	ylabel = "Declination";

	x_options = "ABCGNSTZH";
	y_options = "ABCGNSTZD";

	num_levels = 64;

}

PGPLOT_Renderer::~PGPLOT_Renderer() {

	cout << "PGPLOT_Renderer::~PGPLOT_Renderer Default Destructor" << endl;

}
void PGPLOT_Renderer::set_device(char *dev) {

	device = strdup(dev);

}
void PGPLOT_Renderer::set_numlevels(int nl) {

	num_levels = nl;
}
void PGPLOT_Renderer::set_limits(float low,float hi) {
	low_limit=low;
	hi_limit = hi;
}
int PGPLOT_Renderer::init() {

	
	if (device == NULL) {
		cout << "PGPLOT_Renderer::init using default device " << DEFAULT_DEV << endl;
		set_device((char *) DEFAULT_DEV);

	}

	cpgopen(device);
	ColourMap = new PGPLOT_byr;
	ColourMap->build_table(64);

 	return EXIT_SUCCESS;

}
void PGPLOT_Renderer::set_window(float min_ra_hours, float max_ra_hours, float min_dec_deg, float max_dec_deg) {

	xmin = min_ra_hours * 60.0 * 60.0;
	xmax = max_ra_hours * 60.0 * 60.0;

	ymin = min_dec_deg * 60.0 * 60.0;	
	ymax = max_dec_deg * 60.0 * 60.0;
	draw_frame();

       
}
void PGPLOT_Renderer::set_labels(string t, string xl, string yl) {

	title = t;
	xlabel=xl;
	ylabel=yl;
}	
int PGPLOT_Renderer::draw_frame() {

	cout << "PGPLOT_Renderer::draw() Drawing Member Function" << endl;
	cpgeras();
 	cpgswin(xmin,xmax,ymin,ymax);
	cpgtbox(x_options.c_str(),0.0,0.0,y_options.c_str(),0.0,0.0);
	cpglab(xlabel.c_str(),ylabel.c_str(),title.c_str());

	return EXIT_SUCCESS;
		
}

int PGPLOT_Renderer::build_image(DataContainer *container, unsigned long xdim, unsigned long ydim, vector<float> &ranges) {

	// This is just a Cartesian projection - but we could make this more general by incorporating FITS

	double cdelt1 = (xmax - xmin)/float(xdim);
	double cdelt2 = (ymax - ymin)/float(ydim);
	double hpx_theta = 0;
	double hpx_phi = 0;
	float tr[6],x[2],y[2];
	int rval=EXIT_FAILURE;
	int NC=5;
	
	tr[0] = xmin;
	tr[1] = cdelt1;
	tr[2] = 0.;
	tr[3] = ymin;
	tr[4] = 0.;
	tr[5] = cdelt2;

	vector<float> values(4);
	vector<float> image;
	vector<float> levels(NC);
	vector<char*> labels(NC);
	for (unsigned long j=0; j<ydim; j++) {
		double yval_in_secs = ymin + j*cdelt2 - 0.5*cdelt2;
		for (unsigned long i=0; i<xdim; i++) {
			double xval_in_secs = xmin + i*cdelt1 - 0.5*cdelt1;

			hpx_phi = (((xval_in_secs)/(60.*60.*24.)) * 2*M_PI);
			hpx_theta = (((yval_in_secs)/(60*60*360)) * 2*M_PI);
	//		cout << hpx_theta << " : " << hpx_phi << endl;
			// should check these values
			rval = container->Get_Sky_Val (hpx_phi, hpx_theta, values);
			if (rval == EXIT_SUCCESS) {
	//			cout << i << "," << j << " " << values[0] << endl;
				image.push_back(values[0]);
			}
			else {
				// ???? not sure why this is here
				image.push_back(ranges[0]);
			}
		}
	}
	ColourMap->set_range(ranges[0],ranges[1]);
	cpgimag( (const float *) &image[0],xdim,ydim,1,xdim,1,ydim,ranges[0],ranges[1],tr);
	cpgtbox(x_options.c_str(),0.0,0.0,y_options.c_str(),0.0,0.0);
	for (int cl=0;cl<NC;cl++){
		levels[cl] = ranges[0] + ((ranges[1]-ranges[0])/10.)*cl;
		labels[cl] = new char[6];
		sprintf(labels[cl],"%3.1f",levels[cl]); 
		fprintf(stdout,"%s\n",labels[cl]);
	}
	cpgsci(1);
//	cpgcont( (const float *) &image[0],xdim,ydim,1,xdim,1,ydim,(const float *) &levels[0], NC,tr);
	
//	for (int cl=0;cl<NC;cl++){
//		cpgconl( (const float *) &image[0],xdim,ydim,1,xdim,1,ydim, levels[cl],tr,(const char *) labels[cl],200,400);
//	}
	x[1] = xmin;
	x[0] = xmax;
	y[1] = ymin;
	y[0] = ymax;
	svPlotKnown(x,y,NULL);
	fprintf(stdout,"\n%f,%f", xmin + 600.,ymin+600.);
	cpgwedg("RI",1.0,4.0,ranges[0],ranges[1]," ");
	return EXIT_SUCCESS;
} 
int PGPLOT_Renderer::draw_pixels(DataContainer *container) {

	// First need to get the number of pixels

	int npix = container->ndata;
	
	vector<float> ranges;
	vector<float> values;
	vector<float> vertices;
	vector<float> point;
	vector<float> locus;
	point.resize(2);
	vertices.resize(8);
	
	if (this->hi_limit == -1 && this->low_limit == -1) {
		container->Get_Val_Limits(ranges);
	}
	else {
		ranges.resize(2);
		ranges[0] = low_limit;
		ranges[1] = hi_limit;
	}
	build_image(container,1024,1024,ranges);

	if (npix > 0) {

		cout << "Will attempt to render " << npix << " pixels onto device" << endl;

		for (int pix = 0; pix < npix ; pix++) {

			container->Get_Pix_Val(pix,values);
			container->Get_Pix_Vertices(pix,point,vertices);
			if (verbose == true) {
				cout << "Pixel: " << pix << endl;
				for (unsigned int i=0; i < values.size();i++) {
					cout << " Val: " << values[i] << " Range[" << ranges[i*2] << ":" << ranges[i*2+1] << "]" << endl;
					for (unsigned v = 0; v < 8; v=v+2){
						cout << "[" << vertices[v] << "," << vertices[v+1] << "] ";
					}
					cout << endl;
				}
			}
			/* need to get the colour index */
//			increment = ColourMap->set_range(ranges[0],ranges[1]);
//			colour = floorf((values[0] - ranges[0])/increment);
//			cpgsfs(2);
//			cpgsci(int(colour));
			// draw_pixel(vertices,values[0]);
			// draw_pixel(vertices,values[0]);
			// draw_point(point[0],point[1],1);
			// oplot_healpix_bounds(512,vertices,values[0]); // need to move this to the data container
		}
	}
	else {
		cout << "Container has no pixels (perhaps a function?) " << endl;
	}
	cpgclos();
	return EXIT_FAILURE;
	
}
	
void PGPLOT_Renderer::draw_point(float xin, float yin, int symbol) {
	float x = (M_PI - xin) * 12.0/M_PI * 60. * 60.;	
	float y =  (M_PI/2. - yin) * 180.0/M_PI * 60. * 60.;	
	cpgsci(99);
	cpgpt1(x,y,symbol);
}

void PGPLOT_Renderer::draw_pixel(vector<float>& vertices, float value) {
	
	vector<float> xpts;
	vector<float> ypts;

		/* yes some hard coding here ... should make this much more general but time is pressing */

	int npts = vertices.size()/2;
	xpts.resize(npts);
	ypts.resize(npts);

	for (int i=0;i<npts;i++) {
		xpts[i] = (M_PI - vertices[2*i]) * 12.0/M_PI * 60. * 60.;		
		ypts[i] = (M_PI/2.0 - vertices[2*i+1]) * 180.0/M_PI * 60. * 60.;		
        }
		/* need to calculate the required colour here */

		/* then set the fill area attribute */

		/* then render the pixel */
	cpgbbuf();
	//	cpgsci(1);
	if (verbose == true) {
		cout << "planting polygon: " ;

		for (int i=0; i<npts; i++) {
			cout << "[" << xpts[i] << "," << ypts[i] << "]" << endl;
		}
	}
	cpgpoly(npts,(float *) &xpts[0], (float *) &ypts[0]);
	cpgebuf();
}

void PGPLOT_Renderer::oplot_healpix_bounds(int nside, vector<float>& ranges, float value) {

	vector<double> ramp;
	vector<double> z;
	vector<double> phi;
	vector<double> vv;
	vv.resize(3);

	float minphi = ranges[6];
	float maxphi = ranges[2];
	float minz = cos(ranges[5]);
	float maxz = cos(ranges[1]);
	
 
	int np = 12*nside; // 200 or 2*nside whatever is greater
//        int np=1000;	
	for (int i=0; i<np; i++) {
		ramp.push_back(double(i)/double(np-1));
//		z.push_back(0.0);
//		phi.push_back(0.0);
		
	}

	double halfpi = M_PI/2. ;
	double fn  = double(nside) ; 
	double fn2 = fn*fn ;
	const double con = halfpi / fn / sqrt(3.) ;
	double sj,si,phic,theta;

	for (int regime=0; regime < 2; regime++) { // either equatorial or polar ... star with equatorial
		for (int kk=0; kk < nside ; kk++) {// 0 to nside -1 
			if (regime == 0) { // equatorial regime
				for (int i =0; i<np; i++) {

					// change here to test for range and push onto vectors

//					z[i] = 2./3. - ramp[i] * (4./3.) ; // set up a range of z for 0.666 to -0.666
//					phi[i] = halfpi*(ramp[i] - kk/fn) ; // set up a range of phi
	
					double tempz = 2./3. - ramp[i] * (4./3.) ; // set up a range of z for 0.666 to -0.666
					double tempphi = halfpi*(ramp[i] - kk/fn) ; // set up a range of phi
					if (tempz >= minz && tempz <= maxz) {
						if (tempphi >= minphi && tempphi <= maxphi) {
							z.push_back(tempz);
							phi.push_back(tempphi);
						//	cout << "vert " << z.size() << " tempphi: " << tempphi << " tempz " << tempz << endl;
						//	cout << "vert " << phi.size() << " min phi " << minphi << " maxphi: " << maxphi << " minz: " << minz << " maxz: " << maxz << endl;
						}
					}

				}
			} 
			else {            // polar regime
				for (int i =0; i<np; i++) {
		//			z[i] = (2./3. + ramp[i] * (1.-(kk*kk)/fn2)/3.);  // range of z
		//			phi[i] = kk * con / sqrt(1.-z[i]); // range of phi
					double tempz = (2./3. + ramp[i] * (1.-(kk*kk)/fn2)/3.);  // range of z
					double tempphi = kk * con / sqrt(1.-tempz); // range of phi

					if (tempz >= minz && tempz <= maxz) {
						if (tempphi >= minphi && tempphi <= maxphi) {
							z.push_back(tempz);
							phi.push_back(tempphi);
						//	cout << "vert " << z.size() << " tempphi: " << tempphi << " tempz " << tempz << endl;
						//	cout << "vert " << phi.size() << " min phi " << minphi << " maxphi: " << maxphi << " minz: " << minz << " maxz: " << maxz << endl;
						}
					}


				}
			}
			if (z.size() > 3) {
				for (int quad=0;quad < 4; quad++) {// quadrant
					for (int j=0; j <= regime; j++) { // eqatorial/polar
						sj = pow((-1.),double(j));
						for (int i=0; i<=1 ; i++ ) {
							si = pow((-1.),double(i));
							vector<float> pixel(2*z.size());
							for (unsigned int k=0; k<z.size();k++) {
								zphi2vec (si*z[k], (j+quad) *halfpi + sj*phi[k], vv) ; // converts (z) and (phi) to xyz
								phic = atan2(vv[1],vv[0]) + M_PI; // in [0,2pi]
								theta = asin(vv[2]) + M_PI/2.; // in [0,pi] 
						//		cerr << "vert " << k << " " << "phi: " << phi[k] << " z: " << z[k] << endl;
								// draw_point(phic, theta,1) ; 
								 pixel[2*k] = phic;
								 pixel[2*k+1] = acosf(theta);
								
							}
							draw_pixel(pixel,value);
							pixel.resize(0);
							
						}
					}
				}
				z.resize(0);
				phi.resize(0);
			}
				
		}
	}
}

